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ABSTRACT 

Results of numerical simulations of shock solutions in a geometrical thin accretion disk 
around a Kerr black hole (BH) are presented. Using the smoothed particle hydrody- 
namics (SPH) technique, the influence of the central object is included by means of 
an effective potential, which is known to mimic the external Kerr geometry very well. 
We first present the theory of standing shock formation in accretion disks around a 
Kerr black hole, and show that the results of our numerical simulation agree very well 
with the theoretical results. Using the pseudo-potential, our analysis gives at least a 
semi-quantitative impression into what a true general relativistic computation might 
reveal. The location of the shocks is found to be different for the co-rotating and the 
counter-rotating fiows. Compared to accretion in the vicinity of a non-rotating BH, 
the position of the shock is shifted inwards for prograde accretion and outwards for 
retrograde accretion. We find that the shocks in an inviscid fiow are very stable. We 
also remove the ambiguity prevalent regarding the location and stability of shocks in 
adiabatic fiows. Finally we sketch some of the astrophysical consequences of our findings 

in relation to accretion disks in Active Galactic Nuclei (AGN) and Quasars. 

Key words: Galaxies: active - Black hole physics - Accretion disks - Shock waves - 
Hydrodynamics - Methods: numerical. 
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strong gravitational pull resulting from general relativistic 
effects. Therefore, the specific angular nionientuni of each 
fluid element can remain almost constant and the flow may 
develop steady shock waves for a certain range of parameters 
(Chakrabarti (f990b) and references therein). Shock waves 
change the basic temperature structure of the accretion disk 
and may be responsible for a large number of observed astro- 
physical phenomena of Active Galactic Nuclei (AGN), such 
as soft X-ray production and time variability. 

In the presented paper, we discuss results of numerical 
simulations of shock formation in an accretion disk around a 
Kerr black hole. We use smoothed particle hydrodynamics, 
which has been well tested for one- and two- dimensional 
accretion disks (Chakrabarti & Molteni f993; Molteni, Lan- 
zafame, & Chakrabarti f994), and the pseudo-Newtonian 
potential around a Schwarzschild geometry (Paczyhski & 
Wiita f980). 

We use an effective potential for the central object 
(Chakrabarti & Khanna f992), which mimics the geom- 
etry around a Kerr BH. The Paczyhski- Wiita potential 
mimics the geometry around Schwarzschild Geometry very 
satisfactorily and a considerable body of the recent as- 
trophysical literature, particularly that dealing with ac- 
cretion disks, uses this potential with satisfactory results 
(e.g., Matsumoto et al. (f984), Abramowicz et al. (f988), 
Kato, Honma, & Matsumoto (f993), Lu (f985), Chakrabarti 
(f990b), Chakrabarti (f990a)) A similar potential obtained 
by Chakrabarti and Khanna (f992) is found to reproduce 
the salient features of hydrodynamic flows around a Kerr 
geometry quite accurately. Incorporation of these potentials 
reduces complications in writing numerical codes, without 
sacrificing the physics, in that a code written in a flat ge- 
ometry can be trivially modified to study flow behaviour in 
general relativistic flows. This approach also gives clearer 
insight into the problem, since the contribution from each 
term becomes apparent. We thus believe that our approach 
is well justified and our major conclusions will remain valid 
even when the full general relativistic equations are solved. 
At present, the only SPH solutions known are for spheri- 
cally symmetric Bondi flows (Laguna, Miller, & Zurek f 993). 
The motivation for doing numerical simulations works stems 
from the fact that analytical treatments have been carried 
out only for steady state solutions. Time dependent flows 
need not choose these solutions. Furthermore, the shock- 
solutions found in analytical treatments are not necessary 
stable and may reflect only a transient state. Thus we be- 
lieve that it is essential to study these solutions using fully 
numerical time-dependent codes. 

For physical reasons shocks solutions exist only in a re- 
stricted parameter range in energy and angular momentum 
(Chakrabarti f990b). We show that in the case of co-rotating 
flows, shocks are formed very close to the black hole, whereas 
for counter-rotating flows the shocks form farther away. Re- 
sults of our simulation agree very well with the theoretical 
shock locations which we also present here. 

The plan of our paper is the following. In the next sec- 
tion, we discuss the theory of shock formation in a one di- 
mensional accretion flow around a Kerr BH and discuss the 
regions of parameter space for which shocks are expected to 
form. The analysis presented is time independent . In Section 
3 we describe the results of the time dependent numerical 
simulations. Finally, in Section 4, we discuss the relevance 



of our results in the astrophysical context and present some 
concluding remarks. 



2 THEORETICAL SHOCK LOCATIONS IN 
GEOMETRICALLY THIN ACCRETION 
FLOWS AROUND A KERR BLACK HOLE 

We consider a geometrically thin, adiabatic accretion disk 
flow in the vicinity of a rotating black hole. This analysis fol- 
lows the criteria given by Chakrabarti for steady accretion 
onto Schwarzschild BH (Chakrabarti f990b; Chakrabarti 
f990a). For the accreting matter we use a polytropic equa- 
tion of state: P = Kp'; P and p denote the isotropic pres- 
sure and the matter density, respectively. The adiabatic in- 
dex (7) is assumed to be constant throughout the flow, 
whereas K is related to the specific entropy (s) of the flow. 
We consider the specific entropy and thus K to be constant 
everywhere, except at the shock-locations, where entropy- 
production is allowed. The stationary solution of the prob- 
lem satisfies the conservation of energy, angular momentum, 
and mass, supplied by the transonic conditions at the sonic 
points and the Rankine-Hugoniot conditions at the shock- 
position. Because the radial flow is pressure driven and 
caused by the pull of the BH, not by dissipation, the specific 
angular momentum A remains constant everywhere. Since 
we are interested in the effects of the rotating black hole, 
and in particular the frame dragging, we use strictly geomet- 
rically thin flow of constant transverse height to avoid com- 
plexity. Properties of the shocks are obtained by following 
the general procedures presented in Chakrabarti (f990b). 



2.1 The model equations 

Although we wish to study hydrodynamical flows around a 
rotating black hole, we do not consider it necessary to use 
the exact general relativistic Kerr space-time. Due to the 
complex nature of the geometry around a Kerr black hole 
(which introduces additional problems in writing a general 
relativistic code) we use the pseudo Kerr potential intro- 
duced recently (Chakrabarti & Khanna f992; Khanna & 
Chakrabarti f 992) to obtain the principle properties of the 
accretion flow. In a similar manner to which the Paczyhski- 
Wiita potential (Paczyhski & Wiita f980) mimics a non- 
rotating black hole geometry, the Chakrabarti-Khanna po- 
tential mimics a rotating black hole geometry on the equa- 
torial plane quite accurately, and thus enables us to pursue 
theoretical as well as numerical work using equations in flat 
geometry. The actual deviation of astrophysically relevant 
quantities (such as the specific binding energy) predicted by 
this potential, compared to the exact theory, is very small. 
Therefore, for all practical purposes, our approach is well 
justified. 

In the following we use as the units of mass, velocity 
and distance, the mass of the black hole Mbh, the velocity 
of light c, and the gravitational radius Rg = GMbh I c , re- 
spectively. The radial effective pseudo-potential for the Kerr 
BH takes the form: 
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with the abbreviations 

2ax 



Xo 



±\a\ 



7 
Xl = — 

2 



X* + a^x^ + 2a^x 

Here a denotes the Kerr parameter and x is the non- 
dimensional radial distance from the BH. In this expression 
a positive spin parameter a represents co-rotating accretion 
whereas the negative a together with a negative sign for xo 
represents counter-rotating accretion. The (constant) angu- 
lar momentum A of the accretion flow is always taken to 
be positive. The physical origin of the effective potential 
and the motivation for the several parameters used, are ex- 
plained in Chakrabarti and Khanna (f992). 

With "d the non-dimensional radial velocity and c, the 
sound-speed, the dimensionless energy equation can be writ- 
ten as 



2 ^7- 



+ g{x) 



(1) 



The accretion flow satisfies the equation of mass conserva- 
tion as given by 



m = "dpxho 



(2a) 



where the assumed constant half-thickness is /to/47r. We de- 
fine an entropy sensitive 'accretion rate' M in terms of "d, 
Cs and the polytropic index n = f/(7 — f): 

M = -Dxc^" . (26) 



Note that although the flow is separately of constant entropy 
before and after the shock, entropy increases at the shock 
and so does M ~ rhK"'. 

The physical quantities within the shock region are not 
resolved on a macroscopic spatial scale much larger than the 
mean free path of atoms of the gas. At the shock the follow- 
ing (Rankine-Hugoniot) jump conditions for energy, pressure 
balance and baryon number respectively, are satisfied: 



^4 



£- 



p+ + p+r+ = P- + p-r. 



(3a) 



(36) 



(3c) 



The subscripts '-' and '-|-' refer to the upstream and the 
downstream gas properties; equations (3a-c) give the quan- 
tities after ('-|-') the shock in terms of the known values 
before ('-') the shock. 

An accreting flow, initially subsonic at a large distance, 
has to pass through a critical or sonic point in order to be- 
came supersonic and thus enable shock formation. To select 
the appropriate hydrodynamical solutions the sonic point 
conditions are derived following the general procedure, de- 
scribed by Chakrabarti (f990a): Differentiating the energy 
and mass conservation equations (eqns. f and 2b) and elim- 
inating dcs/dx results in the condition. 



rf»? Cj Cj dg 

dx i} X dx 



(4) 



For a stationary solution with supersonic inner boundary 
condition, the left-hand side and the right-hand side of equ. 
(4) vanish simultaneously at the critical points. The vanish- 
ing of the left-hand side defines a sonic point. 



|»?(Xc)| = Cs(Xc) 



(5) 



Figure 1. Contours of different constant 'accretion rate' AA cor- 
responding to £' = 0.08, A = 1.8, a = 0.95. The transonic solution 
labeled as "SS" represents the physical stationary solution for 
the Ai— = 1.7891 X 10~ , the post-shock solution corresponds to 
A^-l- = 2.5445 X 10~^. The stable shock location is indicated. 



whereas vanishing of the right-hand side gives the sound 
speed at the outer sonic point; the subscript c denotes the 
quantities at the critical points. 

To obtain the full set of solutions which include shock 
waves, the system of equations (1), (2b), (3a-c), (5) must 
be solved simultaneously. We do not repeat the general 
procedure, which has already been discussed in detail in 
Chakrabarti (1990a). Instead, we present a single exam- 
ple of a solution which includes a standing shock. Fig. 1 
presents the contours of constant "accretion rate" M in 
the Mach number (Mach) versus radial distance (log (x)) 
plane. The specific energy is chosen to be £ = 0.08, the 
specific angular momentum A = 1.8. The Kerr parameter is 
a = 0.95. The solution with M- = 1.7891 x 10"* labeled 
as "SS" represents the physical stationary solution which 
passes through the the outer sonic point. The post-shock 
flow with A^-i- = 2.5445 x 10~ passes through the inner 
sonic point. The jump in M from pre-shock to post-shock 
value reflects a jump in the entropy function K. The sta- 
ble shock is located at x = 7.75 in this case (see also figure 
5a). Another shock at x = 2.28 is predicted, but it is eas- 
ily argued that this location is unstable (cp. Chakrabarti & 
Molteni (1993)). 



2.2 The parameter space relevant for shock 
formation 

It is well known that shocks can form if a gas meets an 
obstacle supersonically. In the present case of a black hole 
accretion, the centrifugal barrier acts as an obstacle to the 
transonic flow and shocks form in the accretion disk. How- 
ever, not every set of initial conditions produces a shock. In 
general, the stronger the potential barrier, the higher the en- 
ergy required to "push" the matter through the barrier. For 
the examples presented here, the typical energies required 
are up to one order of magnitude higher for prograde ac- 
cretion than for retrograde accretion (fig. 2). It is therefore 
crucial to have a prior understanding of the regions of pa- 
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Figure 3. Variation of the locations of the shocks with the energy 
as a function of various angular momentum (prograde flows). 



Figure 2. Regions of the parameter space (spanned by (£',A)) 
where stationary shock solutions are predicted. Upper fig. for co- 
rotating, and lower fig. for counter-rotating flows. Black hole ro- 
tation parameter a is as marked. 



raineter space for which stationary shocks may form. Fig. 
2 shows some typical regions in (£, A) space where station- 
ary shocks are both predicted and effectively reproduced by 
our time-dependent simulations. Results are presented both 
for co-rotating and counter-rotating flows and for small and 
large values of the BH rotation parameter a. 

The shocks form closer to the black hole, in the case of 
prograde accretion flows, than they do in retrograde flows. 
Fig. 3 shows the variation of the shocks location with energy, 
for different values of the specific angular momentum. For 
rapidly rotating black holes, the shocks are usually located 
less a few tens of the Schwarzschild radii. On the other hand, 
for slowly rotating black holes (a = 0.4) the shocks are lo- 
cated between a few tens and some hundred Schwarzschild 
radii. Shocks are produced for lower energies and in general 
for higher angular momentum of the gas flow in the case of 
counter-rotating black holes; the stationary shocks are lo- 
cated far away from the black hole horizon (Fig. 4). The 
requirement of higher angular momentum stems from the 



Figure 4. Variation of the locations of the shocks with the energy 
as a function of various angular momentum (retrograde flows). 



fact that some of the centrifugal support is 'spent' to fight 
against the frame-dragging. In general, it is observed that 
for a fixed A the shock strength, given by the ratio of the 
pre-shock to the post-shock Mach number, increases as the 
energy parameter decreases. 



3 RESULTS OF NUMERICAL SIMULATION 
OF ACCRETION DISKS 

We integrate the classical one dimensional hydrodynamical 
equations of a compressible inviscid fluid for axial symmet- 
ric configurations. Let us denote the cylindrical coordinates 
by X (radial distance from the 2-axis of symmetry) and 4>. 
The assumed axial symmetry demands ■^ = 0- With the 
Lagrangian derivative 



+ v -y 



D_ _ d_ 
Dt ~ dt 

we have for the mass conservation 
D Q Q d{xva 

Dt X dx 
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For the radial and the aziniuthal nionientuni we use 



Dt 



ldp_ 
Q dx 



dgjx) 
dx 



In 



V^Vx 



where g{x) describes the effective potentiai of the Kerr BH 
as given in sec. 2.1. The energy equation which describes the 
behavior of the internal energy per unit mass e is 



Dt 



1 d{xvx 
X dx 



The procedure for iniplenientation of the above equations 
into the SPH fornialisni is presented in Chakrabarti & 
Molteni (1993) and is not repeated here. 

The general procedure followed in our simulations is 
presented in Chakrabarti and Molteni (1993); there they 
showed that the time dependent code reproduces — exactly, 
within the numerical resolution — the analytical solutions. 
We use a similar SPH code except that the effective pseudo- 
Kerr potential (Chakrabarti & Khanna 1992) is used instead 
of the pseudo-Schwarzschild potential (Paczyhski & Wiita 
1980). The boundary conditions are as follow: At the in- 
ner radius (r = 1.5 x GMth/c ) the smoothed particles are 
simply cut away from the computation; since the flow here 
is extremely supersonic no feedback disturbances appear on 
the upstream flow. The smoothed particles are introduced 
with the azimuthal velocity v,f, = X/x, at an injection po- 
sition (xi) downstream of the outer sonic point where the 
radial gas speed is already supersonic. The injection velocity 
Vj;{xi) and the sound speed Cs{xi) are chosen at this outer 
boundary according to the prescription for which the theo- 
retical steady solution is known from the theory presented 
in sect. 2.1. A new particle is injected every time an injec- 
tion sphere S (around the injector position) is free. The flow 
is smooth for small injection sphere (since many particles 
overlap). For greater injection areas the flow is more turbu- 
lent. We note that for extremely smooth flows the solution 
obtained in these simulations remains supersonic through- 
out according to un-shocked solution. When the flow at the 
outer boundary is more turbulent, the flow switches to the 
solution which includes a shock. In this case, the flow passes 
through both the sonic points. 

Figures 5, 6 show the shock profiles for different param- 
eters a, £, A. To within the range of resolution 2h (h being 
the particle size used in the simulation) the shock positions 
and strengths agree very well with the theoretical results. 
All the temperature profiles are given in nondimensional 
units of the thermal energy per unit mass. To estimate the 
temperature values in Kelvin for this one-dimensional com- 
putation an appropriate accretion rate must be chosen and 
the vertical extent of the flow must be assumed. The accre- 
tion rate is expressed in units of the Eddington accretion 
rate Me ~ 0.2MsMQyear~ , where Ms denotes the mass 
of the BH in units of 10 solar masses (Blandford 1990). If 
we assume the flow at the inflow region to be in vertical 
force equilibrium and use the resulting vertical extent for 
the whole radial region, we obtain a temperature presented 
in figs. 9. If the radiative pressure Pr is much higher than 
the gas pressure Pg we can assume that the temperature 
of the ID disc is given by T^ = P3c/4:Cr, where the total 
pressure P = Pr + Pg is given by the polytropic pressure 
P = (j — l)eUth- Here the mass density q = pQo is com- 
puted from the nondimensional density p (which is obtained 



Figure 7. Theoretical and numerical Mach and temperature 
profiles. Temperature profiles are scaled in units of the ther- 
mal energy per unit mass. The parameters are a = 0.95, £' = 
0.08, A = 1.8 with the theoretical shock position close to the BH 
at Xs = 7.75. The theoretical shock position within the range of 
resolution is marked by the two vertical lines. The right parts 
of both figures show some enlargement of the shock-region. The 
corresponding temperatures in Kelvin, dependent from the mass 
accretion (-M) in units of the Eddington accretion rate and the 
mass of the BH (Mj^jj) are displayed in fig. 9, left. 



from the computation) and scales with the properties of the 
accretion flow at the injection point (xi): 



go 



M 



4:TrXtZdVj:{Xt 



with Zd being the assumed vertical disk equilibrium extent 
at Xj. Uth is the thermal energy per unit mass which results 
from the analytical as well as from the numerical solution, 
7 = I denotes the adiabatic index and cr denotes the Stefan- 
Boltzmann constant. 

Figures 7, 8 show the excellent agreement between the- 
oretical and numerical Mach number and temperature pro- 
files. These two cases have the parameters (fig. 7) a = 
0.95, £ = 0.08, A = 1.8 (see also fig. 1), and fig. 8 a = 0.4, £ = 
0.008, A = 3.6. 

Figure 9 shows the corresponding temperature in Kelvin 
for the choice of various black holes masses and accretion 
rates. The shock positions agree to within the spatial reso- 
lution of the particle size h. The temperature and the Mach 
number agree exceptionally well in the subsonic branch; only 
in the supersonic branch is there a small difference that 
it is clearly due to the visible jump at the inflow bound- 
ary. The shock starts forming near to the centrifugal barrier 
and travels backwards with time and eventually become sta- 
ble practically at the theoretically predicted shock location 
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Figure 5. Prograde flows: Variation of Mach number, temperature profiles (scaled in units of the thermal energy per unit mass), 
and non-dimensional mass density for different input parameters a, £', A. The shock positions and strengths agree within the range of 
resolution (2/i) with the theoretical predictions. In each figure theoretical shock position with the range of resolution i^Xg i 2h) is marked 
by the two vertical lines. 



Steady state shocks in accretion disks around a Kerr black hoi 



Figure 6. Retrograde flows: Variation of Mach number, temperature profiles (scaled in units of the thermal energy per unit mass), 
and non-dimensional mass density for different input parameters a, £, X. The shock positions and strengths agree within the range of 
resolution {2h) with the theoretical predictions. In each figure theoretical shock position with the range of resolution (xg i 2h) is marked 
by the two vertical lines. 
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Figure 8. Theoretical and numerical Mach and temperature pro- 
files. Temperature profiles are scaled in units of the thermal en- 
ergy per unit mass. The parameters are a = 0.4, £' = 0.008, A = 
3.6 (predicted shock position far away from the BH horizon at 
Xs = 80.66). The theoretical shock position within the range of 
resolution is marked by the two vertical lines. The right parts 
of both figures show some enlargement of the shock-region. The 
corresponding temperatures in Kelvin, dependent from the mass 
accretion (-M) in units of the Eddington accretion rate and the 
mass of the BH {M^h) ^re displayed in fig. 9, right. 



(c.f. Chakrabarti & Molteni (1993)). In the simulation, we 
chose a = 1, /? = 2 for the artificial viscosity parameters 
to smooth out the oscillations at both sides of the shock 
(though a little overshooting is still visible). With a higher 
artificial viscosity, the shock seems to become diffusive and 
spreads out on both sides. 

In passing, we remark that we have performed several 
simulations for which theoretical treatments suggest that no 
shock forms. We do not find any shock formation in the nu- 
merical simulations either. Even when large perturbations 
are introduced, the flow seems to return to its steady su- 
personic solution after a short period of transients. One of 
the most significant results of our simulation is that the flow 
sfearfy becomes asymptotically and smoothly passes through 
hoth the saddle type sonic points simultaneously, without 
ever shuttling back and forth and causing bi-periodicity as 
one may naively expect . 



4 OBSERVATIONAL SIGNIFICANCE AND 
CONCLUDING REMARKS 

In this paper, we have studied the nature of the Rankine- 
Hugoniot shocks in thin, axisymmetric accretion around a 
rotating black hole both theoretically and through numer- 



ical simulation. We find the agreement between these two 
approaches to be excellent. This provides us with confidence 
that the code can be trusted not only for thin flows, but for 
more general accretion flows as well. 

Our numerical simulations indicate that shocks may be 
common in accretion disks, especially if the angular momen- 
tum of the disk at a large distance is low compared to the 
'Keplerian' value: the angular momentum of the circular or- 
bit around a Kerr BH. This is a possibility if the accretion 
takes place from winds in star clusters. We have also shown, 
that out of two theoretically predicted shock locations, only 
one is stable, as Chakrabarti & Molteni (1993) showed for a 
nonrotating BH. This removes the ambiguity in shock loca- 
tions in adiabatic flows. 

It is generally believed that the black holes in galactic 
centers are rotating. One of the ways to be certain that a 
black hole is rotating is to look for features in a disk which 
are sensitive to the spin of the black hole. As Chakrabarti 
(1990a) has pointed out, a stable standing shock could be 
one such feature. The shocks in co-rotating disks are closer 
to the hole, which means that the post-shock region of the 
disk is not very big. In this case, perturbations in the ac- 
cretion rate should induce faster variabilities in the spectra 
from the disks, than for accretion by a Schwarzschild BH. 
Furthermore, the location of the shock itself is very sensi- 
tive to the energy and angular momentum of the flow. We 
suspect that these variabilities may cause the X-ray flicker- 
ing that is observed in AGNs. On the contrary, the shocks 
in counter-rotating disks form very far away from the BH, 
causing a significant part of the disk to be hotter than a Ke- 
plerian disk, producing excess UV and possibly soft X-rays. 
The variability time scale should be longer in this case. 

Counter-rotating flows are possible in the case of wind 
accretions (re-accretion of earlier expelled low-angular mo- 
mentum winds). A particularly interesting system to observe 
would be one, where the variability time scale switches (at 
least once!) from a short time scale to a long time scale 
(or vice versa). This could represent a case of wind accre- 
tion switching from a predominantly co-rotating manner to 
a predominantly counter-rotating manner (or vice versa). 
It may be quite possible to observe this, since even a small 
change in the direction of the wind at a large distance might 
change the character of the accretion from co-rotating to 
counter-rotating. 

Shocks in accretion disks may be useful in accelerat- 
ing particles (by, say, first order Fermi processes) and thus 
producing cosmic rays. It is possible that details of the 
resulting cosmic ray spectrum strongly depends upon the 
Kerr parameter a. If scatterers (e.g. weak magnetic field) 
are present close to the shock front in the pre-shock and 
the post-shock region, particles cross the shock front several 
times and gain momentum from the scatterers. As the result- 
ing energy spectrum is mainly determined by the difference 
between the upstream and downstream velocity of the scat- 
tering centres, the momentum (p) distribution function of 
the scattered particles is roughly f(p) oc p~ '"'('"" > where r 
is the shock compression ratio, the ratio of the upstream to 
the downstream velocities (Blandford 1990). 

We also note that in the case of retrograde accretion, 
the spatial efficiency of Fermi acceleration of cosmic rays in 
the standing shock could be higher compared to the same 
process in a co-rotating BH, or even a Schwarzschild BH, 
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Figure 9. Temperature profiles in Kelvin for the simulations presented in the figs, different masses (Mj^jj) of the BH and different mass 
accretion rates (M) related to the Eddington accretion rate Mjt; (cp. sect. 3). For both figs, the parameters are for the graphs from up 
to down: Mbh = 6Mq: M = O.IMe, M = O.OIMe and Mbh = ICMq: M = O.IMe, M = O.OIM^. 

just due to the fact that the shock fronts are more extended This paper has been produced using the Blackwell Scientific 

in these cases. The detailed dependence of the strength of Publications I^TjtX style file, 

the shocks on the BH parameters especially for more real- 
istic 2D-simulations, and details of the emitted cosmic ray 
energy spectrum, will be a topic for future investigation. 
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